COUPLED MODE EQUATION MODELING FOR OUT-OF-PLANE GAP SOLITONS IN 

2D PHOTONIC CRYSTALS 



TOMAS DOHNAL* AND WILLY DORFLER* 

Abstract. Out-of-planc gap solitons in 2D photonic crystals are optical beams localized in the plane of periodicity of the 
medium and delocalized in the orthogonal direction, in which they propagate with a nonzero velocity. We study such gap solitons 
as described by the Kerr nonlinear Maxwell system. Using a model of the nonlinear polarization, which does not generate higher 
harmonics, we obtain a closed curl-curl problem for the fundamental harmonic of the gap soliton. For gap solitons with frequencies 
inside spectral gaps and in an asymptotic vicinity of a gap edge we use a slowly varying envelope approximation based on the 
linear Bloch waves at the edge and slowly varying envelopes. We carry out a systematic derivation of the coupled mode equations 
, (CMEs) which govern the envelopes. This derivation needs to be carried out in Bloch variables. The CMEs are a system of coupled 

■ nonlinear stationary Schrodinger equations with an additional cross derivative term. Examples of gap soliton approximations are 

fj^ ' numerically computed for a photonic crystal with a hexagonal periodicity cell and an annulus material structure in the cell. 
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1. Introduction. Maxwell's equations for electromagnetic waves in Kerr nonlinear dielectric materials 

read 

a t D = Vx H, (1.1a) 
fi d t H = -V x £ , (1.1b) 

rn! v-x> = o, (i.ic) 

oo ' 
in 
m 



V-H = (l.ld) 

for the electric field £, magnetic field "H, the electric displacement field T> with the constitutive relations 

V = e (n 2 £ + V NL ) , 

3 ^ 

^NL,i = X?jlm £ & £ ™ fori = 1,2, 3. 



j,Z,m— 1 



£0) Mo are the electric permittivity and magnetic permeability of vacuum, respectively, x h4 nix) is the refractive 
index of the medium, and x H> x^\ x ) ls t ne cubic electric susceptibility of the medium. 

We consider a 2D photonic crystal, i.e. we assume that the material coefficients change periodically 
on a plane and are independent of the orthogonal component on that plane. Let be linearly 

independent lattice vectors defining the Bravais lattice A := span z {a'^, a^} of the crystal. Then the required 
periodicity reads 

nix) = nix + R) e R, 

(13) 

x {3) (x) = X ( 3 \ X + R) e R 3x3x3x3 for all x G M 3 and all R E A. 

Without loss of generality we assume that the crystal is homogeneous in the a^-direction, i.e. = = 
and d X3 n — d X3 Xiji m = ® f° r an hjih m - We denote by U the Wigner-Seitz cell corresponding to the Bravais 
lattice. We use b^\b^ to denote the pair of vectors satisfying a™ • b^> = 2TrSi,j for i,j E {1, 2}, and let the 
reciprocal lattice be A* := span^jM 1 ), b^}. B denotes the first Brillouin zone, i.e. the Wigner-Seitz cell of 
the reciprocal lattice. 
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Note that from the relations in (|1.2j) it is clear that we are neglecting losses, material dispersion as 
well as higher order nonlinearities and assuming that the third order nonlinear response of the medium is 
instantaneous. 

We will consider monochromatic waves propagating in the ^-direction, i.e. waves propagating out of the 
plane of periodicity of the 2D crystal, and use the ansatz 

(£, H, V)(x, t) = eK^a-wt) (£ ; h 7 D)(x u x2\u) + cc, (1.4) 

where k£I and cc. denotes the complex conjugate of the first term on the right. The ansatz (|1.4j) contains 
no higher harmonics, which is valid if the above form of Pnl is replaced by a time averaged one, see below. 
Alternatively, a physical justification of neglecting higher harmonics is based on the lack of phase matching 
and absorption. 

Note that for the field (|1 .4[) the divergence free conditions (jl.lc[) and (jl.ld[) are automatically satisfied 
provided w/0 since the spatially dependent parts 

(£,%T>)(x;u) := e iKX3 (E, H, D) (s x , x 2 ; u) 

satisfy 

V = -V x H and = --V x £, (1.5) 

LO UJ 

and thus V • T> = V • % = 0. Since our analysis below is for gap solitons with uj close to a band edge, the 
condition ui ^ is for us restrictive only when uj — is in a gap and lies near a band edge. Note also that 
even if higher harmonics are accounted for, the divergence free conditions are still satisfied for u =/= as (|1.5[) 
then holds for each generated harmonic. Clearly, only odd, i.e., (2n + l)-th, n G Z, harmonics are generated. 
We will assume a centrosymmetric and isotropic x^-tensor, which leads to the simplification 

Pnl =X<f 

where Xd? '■— Xmi = X2222 = X3333 f° r : ( x ii x 2) £ K 2 -> 1, see [2TJ Sec. 2d]. Inserting the ansatz (|1.4p 
in the nonlinearity Pnl clearly generates the harmonics e ± 3l ( Ka:; 3-wt)_ These are, however, typically neglected 
based on the physical arguments that the fundamental harmonics e ±i( rex 3 — w *) and the higher harmonics are 
not phase matched and that at the higher values of frequency (i.e. at ±3w) material absorption is usually 
large preventing the generation of significant fields at these frequencies, see e.g. [9]. Considering only the 
fundamental harmonics, the nonlinear polarization for the ansatz (|1.4|) becomes 

Pnl = xf {^\E\ 2 E + E ■ EE)e i{KX3 - ut) + cc. , (1.6) 

i.e. 

rq s /(3|Bi| 2 +2|B2| 2 +2|E 3 | 2 )£;i + (£; 2 2 + E3 2 )Ei \ 

Pnl - Xci (2|B 1 | 2 +3|£ 2 | 2 +2|£ 3 | 2 )£ 2 +(B 2 +B 2 ) J e 2 e'( ras - ui '+c.c.. (1.7) 

\ {2\E 1 \ 2 +2\E 2 \ 2 +3\E 3 \ 2 )E 3 + {El + El)E 3 J 

Another widely used model for the nonlinear polarization is 

where [/] av denotes the time average of / over the period of /, i.e. over t € [0, 7r/w] for / = £ ■ £ , cf. [2l)ll2"7] . 
The averaging generates no higher harmonics so that in this model (| 1 . 6|) is exact. Note that the Kerr nonlinear 
problem including all higher harmonics has been recently considered for a ID periodic structure in [25] . 
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In the following we rescale the frequency by defining 

~ w 
to := — 
c 

but drop the tilde again for better readability. For convenience we will denote the square of the refractive 
index by 

rj(x) := n 2 (x) for all x E M 3 . 
With the ansatz ()1.4[) equations (|l.la|l and (|l.lb|) become 

-icwD = Vxff + i(oj x H, 
icuj^qH = V X E + i f o J x E. 

Since all our functions are independent of £3, we let from now oni= (^1,^2) G Using the fact that E 
depends only on x\ and X2, a second order formulation of (|1.8[) reads 

(L-lo 2 t 1 )E = uj 2 P nl , (1.9) 

where 

LE :=VxVx£ + k( d7 2 E 3 \ + n 2 ( eI) , (1.10) 

\d xl E 1 +d X2 E 2 J ^ / 

and 

^L=x£°(2|£| 2 £ + 

Having determined the magnetic field can be recovered by 

H=-^(vxE + i(l)xE). 

Based on the analogy with the periodic nonlinear Schrodinger equation [53], equation (ll.9[) is expected 
to have localized iJ(curl, Resolutions E for any w in a spectral gap of the linear problem Lu = u> 2 iju. Such 
solutions are called gap solitons. The aim of this paper is to provide an approximation of gap solitons E of 
(|1.9[) for lu in an e 2 - vicinity (0 < e <C 1) of a gap edge using a slowly varying envelope approximation. As 
we show, envelopes of such gap solitons satisfy a system of nonlinear constant coefficient equations, so called 
coupled mode equations (CMEs) posed in the slow variables y = ex. The CMEs can be numerically solved 
with less effort than the nonlinear Maxwell system (|1.9p in the variable x. An asymptotic approximation of 
a gap soliton of (|1.9[) near a gap edge is then the sum of linear Bloch waves at the edge, modulated by the 
corresponding envelopes. 

Asymptotic approximations via CMEs have been analyzed for gap solitons of the stationary periodic 
nonlinear Schrodinger equation in ID [51] as well as in 2D [T5] [TJ]. In these works the approximation 
via CMEs was also rigorously justified using Lyapunov-Schmidt reductions. Gap solitons of the nonlinear 
Maxwell's equations have been approximated by CMEs in the case of ID photonic crystals with a small 
(infinitesimal) contrast in the periodicity [THIHHHS], where [TB] considers gap solitons modulated also in time. 
To our knowledge the problem of a systematic CME approximation of gap solitons of nonlinear Maxwell's 
equations describing 2D or 3D photonic crystals does not appear in the literature. Although CMEs have been 
formally derived for pulses in Maxwell's equations with a 2D periodic medium of small contrast [2][T][Tl], these 
pulses cannot be true gap solitons because in 2D and 3D a large enough contrast is necessary for the opening 



(1.8a) 
(1.8b) 
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of spectral gaps. In this paper we consider a 2D photonic crystal with a finite contrast in the periodicity. For 
our examples we use a photonic crystal which has several spectral gaps [3]. 

Besides the above cited works on coupled mode modeling of gap solitons there are a number of papers on 
the slowly varying envelope approximation of nonlinear pulses in periodic structures with the pulse frequency 
lying within the spectral bands. The envelope in this case can be typically modeled by the time dependent 
nonlinear Schrodinger equation and the approximation holds on large but finite time intervals [HI IE1 HH] • 

The rest of the paper is organized as follows. In Section [2] we study the linear band structure w n (fc) of 
(jl.9[) (with x c i = 0) and obtain thus the linear spectrum of the problem. We also discuss possible symmetries 
in the band structure and among the corresponding Bloch waves. An example of a photonic crystal from [3] is 
then provided, for which the band structure is numerically computed and three band gaps are observed on the 
positive half axis uj > 0. In Section [3] we present a slowly varying envelope approximation of gap solitons of 
(ll.9[) for cj in the vicinity of a spectral edge and carry out a systematic formal derivation of CMEs describing 
the envelopes. Next, examples of CMEs are presented for the concrete photonic crystal given in Section [2] 
as well as for other theoretical situations. Here the symmetries in the band structure and among the Bloch 
waves play an important role in determining properties of the CME coefficients. In Section UJ we plot the 
approximation of two gap solitons in the chosen photonic crystal. The approximation requires computing the 
Bloch waves at the edge and solving the corresponding CMEs. 

2. Linear Band Structure. 

2.1. The periodic eigenvalue problem. We study first the linear problem 

Lu = u 2 r)u on K 2 (2.1) 

and define the band structure as well as the linear Bloch waves. 

By the Bloch-Floquet theory, see [JJ5] or [T5J Ch. 3] , solution modes of (|2.1[) are given by the Bloch waves 
u n (k; . ) for n £ N that satisfy 

Lu n (k; . ) = u> n (k) 2 r)u n (k; . ), 
u n (k; . + R) = u n {k; .)e ihR for all R e A, 

where k = {ki, ^2) sweeps the first Brillouin zone BcR 2 . 

It is well-known that L is self-adjoint and has a compact inverse and that there thus exists a sequence of 
eigenvalues {w n } n >i with limn-^ uj n = 00 and each eigenspace is of finite dimension. These eigenvalues are 
nonnegative and we use the natural ordering w„-i < uj n for n > 1. The mapping k 1— > uj n (k) is called the n-th 
band of the spectral problem (|2.2[) . Of course, (|2 .2[) allows also non-positive bands —u> n . These are typically 
labeled via w_„ = -w„ and will play no role in our analysis. We therefore restrict ourselves to uj n > for 
n£N. The Bloch waves in (|2.2p can be written in the form 

u n (k;x) = p n (k;x)e lk ' x , 

where the p n are A-periodic in x, i.e. p n {k; x + R) = p n (k; x) for all x 6 U, R £ A. These satisfy the eigenvalue 
problem 

L(k) - ujl(k)ri(x) \ p n (k; x) = for all x £ U, 

(2.3) 



with 



p n (k; x + R) = p n (k; x) for all x £ dU and all R £ A, 



L(k)p n (k] x) = (V + ik') x (V + ifc') x p n (k- x), 
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where k = (fci, £2) S B, k' = (k%, k2, n) T . Since p n is a^-independent, L(k) can be written as 

/ K 2 -(Sx 2 +ifc 2 ) 2 (dx, 1 +iki)(d X2 +ik 2 ) i«(9 a;i +ife 1 ) \ 

L(fc) = (9» 1 +ifei)(9 ! » 2 +ifc2) K 2 -(9 xl +ife!) 2 i K (9 X2 +ifc 2 ) 

\ ^(ao-j+ifex) iK(9 X2 +i/c 2 ) -(a xl +ifei) 2 -(a X2 +ifc 2 ) 2 / 
In the variable k the Bloch waves u n and the eigenvalues w„ are easily proved to fulfill 

uj n (k) =uj n {k + K), p n (k + K;x)=p n (k;x)e- iK x for all x G U, K G A*. (2.4) 
Due to the self-adjoint nature of L(k) we can normalize the Bloch functions via 

(p„(k; . ), r)p m (k; . )) = 6„ jm , (2.5) 

where (f,g) = (f,g) L ^uy 3 = lu f( x ) ■9( x ) dx for f>9 ■ K 2 -> C 3 . 

For purposes of the later asymptotic analysis of gap solitons we also present calculations of first and second 
order derivatives of the bands at extremal points. Suppose the band w n< , has an extremum at k = fc* G B and 
denote :— w n<- (fc*). By direct differentiation of ()2.3[) we see that the "generalized Bloch functions" d kj p n ,, 
for j G {1,2}, are solutions of the system 

(Z(K) -uj^i]j d k:l Pn, (k*; •) = -d kj L(k*)p nr (k*; .). (2.6) 

Applying the differentiation d ki kj , for i, j G {1, 2}, to (|2 .3p and evaluation at n = ra*, k — fc* yields 

-^*»7(a;)) dl z kj p n ,(h;x) 

- d kt L(h)d kj p ru (K; x) ~ d kj L(k*)d ki p nt {K;x). 

Necessarily, due to the Fredholm alternative, the right hand side is L 2 -orthogonal to p n ,(k*; ■ ), which lies in 
the kernel of L(fc*) — oj^t] with periodic boundary conditions on U. This yields the formula 

(0fcWn.(fc*))ij = d k z ,k 3 ^nA k *) 

= Tj^j- (dl i k] L(k*)p n , (k*\ . ) + d k jL{k*)d k .p n ,{K\ . ) + d kj L(k«)d kl p„,(k*; . ),p n ,{k*, -)^ ■ (2.7) 
A straightforward differentiation of L(k) yields 

„ / i(dx 2 +ifc*,2) K 

d kl L{K) = i(a X2 +ifc», 2 ) -21(8^+1^,1) 

\ — K -2i(9 xl +ifc,,i) 

/ -2i(S X2 +ifc„. 2 ) i(9 xl +«c*,i) 

d k2 L(h) = i{d xl +ik.,!) 

\ — k -2i(d X2 +ifc„, 2 ) 

where for j G {1, 2, 3}, is the j-th component of fc*. With these the explicit forms of (|2.7[) read 

I If i(dx 2 +ifc*,2)9fc 1 Pn„,2(*:»; • )-t<? fcl p„ t >3 (fe* ; . ) 

f i/s*,2)9fc 1 Pn»,l(^»; • ) — 2 K d a>l+ i k*,l)dl°lP''*,2{ fl> *; ■ ) + 

-2i(a xl +ifc»_ 1 )ai. 1 p„, i3 (fc»; . )-KO fcl p„»,i(fc t : . )+p„ t>3 (fc„; . ) 



I I / i(3 X2 +ifc. i 2 )9 fcl p„„ . 2 (/c»; . )-Ka fcl p„ t>3 (fc*; . ) \ \ 

u n , (K) = — ( i(e X2 +ife ti2 )a fcl p„,, 1 (fe»;.)-2i(a xl +ifc», 1 )a t;i p„,, 2 (fc»;.)+P„»,2(fc,;.) ) ,p n ,{K; .)), (2.9) 
\ V -2i(^ 1 +ifc», 1 )9 fe ,p n ,, 3 (fe»;.)-«9jfc 1 p»»,i(S!,;.)+p n «,3(fc.;-) / / 



^ / / -2i(a x2 +ife». 2 )9 fc2 p„»,iC=»; • )+i(<?x 1 +»c,,i)9fc 2 P>«»,2(fc,; ■ )+p„..i (fc. ; ■ ) ' 

dluj n ,(h) = — ( i(9» 1 +ife,,i)8* 2 Pn„i(fc*;-)-«9* 2 fn», 3 (fc.;.) ),Pn„0*;-))> ( 2 - 10 ) 

w * \V -K9fc 2 p n »,2(fc.;.)-2i(a x2 +ifc,, 2 )a fc2 p„,, 3 (fc»;.)+P, lt ,3(fc.;.) 
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and 



/ / -2i(a x2 +ife» i2 )9fc 1 p„»,i(fc,; . )+i(d xt +it,,i)4 1 p„, , 2 (fc»; ■ )+i(<3 X2 +ik»,2)3fc 2 P™» ,2(^*1 • ) 
d 2 k W n , (fc*)= 75^— ( ( i (^i+ ifc »,i)9fc 1 p n ,,i( fc »; ■ )+i(9:c 2 +ife»,2)9fc 2 Pr 1 ,,i(fc*; ■ )- 2i(a x j , 1 )d fc2 p„„ ,2 (*:, ; • ) 
' 2 * \\ -2i(9x 2 +ifc* 1 2)afc 1 P, 1 »,3(fc.;-)-2i(a xi +ifc»,i)a fc2 p„,, 3 (fc,;.) 



-2i(9„. 2 +ifc» i2 )9 fc p„, i3 (fe»; . )-2i(d xl +ife» 1 i)9| S . 2 p„„ , 3 (fc» ; . ) 

-K3fc 2 p„„,3(fc»; ■ )— Pn,,2(fc»; ■ ) \ \ 

-tS* 1 P»»,3(fc.;-)-p».,i(*.;-) ) ,p nt (k*; .) 

— «(9ls 1 Pn,,2(fc»; • )+9* 2 P«»,l(' s *i ■ )) / / 



(2.11) 



2.2. Symmetries of the Band Structure and the Bloch waves. Symmetries in the refractive index 
function rj yield symmetries in the band structure and among Bloch waves. We restrict our attention to the 
cases of discrete rotational and axial reflection symmetry, which are relevant for the example we present below. 
The results of this section will be important when determining properties of the coefficients of coupled mode 
equations in Section [3.41 

2.2.1. Rotational symmetry. Assume that the photonic crystal satisfies the rotational symmetry 

T)(x) =r)(r a (x)) for all x G R 2 (2.12) 
for some a G (— 7r,7r] with the rotation r a defined by 

r (x) = ( cos ( a ) x i- sin ( a ) x 2 \ 

a V / \ sin(a)xi+cos(a)a:2 J 

Below we use the notation r a (v) — (cos(a)ui — sm(a)v2, sin(a)t>i + cos(a)«2) T if w is a two dimensional vector 
v <E C 2 and r a (v) — (cos(a)ui — sin(o;)u2 ) sin(a)ui + cos(a)^2, «3) T if v is a three dimensional vector v G C 3 
The symmetry (|2.12l) implies a symmetry of the Rayleigh quotient corresponding to the eigenvalue problc 
} and thus a symmetry of the band structure. In detail, for k G B we have 

. J^KV + ifcQ x w(x)\ 2 dx 

0J n (k) = mm max — — , 

V<ZH™\U) «>eV>/o Jjj^xjlw^x)]- 1 dx 

dim V—n 



3 



tion 



and the corresponding extremal point is p„(fc; . ). Due to the relati 

((V + ir a (k')) x /) (r a (x)) = r a [(V + ifc') x r_ Q {f(r a (x)))} for all smooth / : 

we get 

/ |(V + ir a (k')) x w(x)\ 2 dx = [ |(V + ife') x r_ Q (w(r a (x))) \ 2 
Ju Ju 

and symmetry (|2.12l) yields 

ri(x)\w{x)\ 2 dx = I rj(x)\r- a (w(r a (x)))\ 2 dx 
*U Ju 

As a result we obtain that 



D 2 . Tn> 3 



da;. 



u n (k) = u n (r a (k)) for all n G N and all k G B. (2.13) 

If Lo n (k) has geometric multiplicity one as an eigenvalue of (|2.3p . we have also a symmetry of the corresponding 
Bloch functions, namely 

Pn{r a (k); x) = e la r- a (p n (k;r a (x))) for all n G N and some a = a(n) G K. (2-14) 
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Note that a renormalization of p n {f a (k),x), in order to obtain a — in (|2.14[) . is in general impossible when 
i" a {k) = k, where k = I reads u k congruent to I" and means k — l + K for some K £ A*. This is because in this 
case p n (Ta(k); x) and p n (k; x) are related by (|2.4[) and a renormalization of the left hand side of (12.14[) would 
affect the right hand side in the same way. When r a (k) is not congruent to k, e.g. when k £ int(B) \ {0}, then 
one can set a — in (|2.14jl . 

^From the symmetry 13|) we can deduce a symmetry of the second derivatives of to n . Using the identity 
dkOJ n (k) = dk(uj n (r a (k))) — (r Q ) T (9feW„)(r Q (fc)), we get by further differentiation 



' 9^w„(r a (A;)) \ j cos 2 (a) sin 2 (a) 
dl 2 uj n {r a (k)) 

\9k u k 2 u n{ra(k)) J 



- sin 



(2a)\ / d 2 k u n (k) \ 



sin 2 (a) cos 2 (a) sin(2a) d k2 w n {k) 
sin(2a) — A sin(2a) cos(2a) J \pl i k2 uj n (k) j 

for all k e B and neN. 

2.2.2. Reflection symmetry. If the photonic crystal satisfies the reflection symmetry 

rj(x) = ri(Si(x)) for all x e M 2 , where Si(x) = (-Xi,x 2 ) T , 

then similarly to Section ?2. 2 .11 we have 

w n (fc) = uj n (~ ki, k 2 ) for all k G B and n £ N. 

Again, if uj n (k) has geometric multiplicity one as an eigenvalue of (|2.3p . then 

p n (Si(k); x) = e la Si (p„ (k; Si(x))) for all n £ N and some a = a(n) £ M. 



(2.15) 



(2.16) 
(2.17) 
(2.18) 



where Si(v) — (—vi,V2,V3) T for v £ 
symmetry (|2.17[) implies 



3 . Just as above, unless k = S\ (k), we can set a = in (|2.18[) . The 

(2.19) 



d 2 kl uj n {k) = (dl 1 uj n )(-k 1 ,k 2 ), d 2 k2 u: n (k) = (dl 2 oj n )(~k 1 ,k 2 ) 
d 2 klM u n (k) = -(dt lM u n ){-k u k 2 ) 



for all k £ B and n £ N. 

An analogous discussion, of course, applies for the reflection symmetry r/(x) = ^(^(a;)) for all x £ R 2 , 
where Szfa) — (xi, —x 2 ) T - One the obtains 



d 2 kl u n {k) = (d 2 kl u n )(k l7 -k 2 ), d 2 k2 uj n {k) = (^ k2 u; n )(h, -k,), 
dfc 1)fe2 w„(fc) = -{dl lM u n ){k u -k 2 ) 
for all k £ B and n £ N and if ui n (k) has geometric multiplicity one as an eigenvalue of (|2.3|) . then 
Pn(S 2 (k)- x) = e la S 2 (p n (k; S 2 (x))) for all n £ N and some a = a(n) £ K. 



(2.20) 



(2.21) 



2.2.3. Combination of rotational and reflection symmetries. If both the reflection symmetry 
(|2.16|) and the rotational symmetry (|2.12[) for some a £ (— tt, — tt], |a| ^ 7r/2, hold, then for k along the 
rays with angles tt/2 — a/2 and — (tt/2 + a/2) the mixed derivative d ki k2 uj n (k) can be expressed in terms of 
d ki uj n {k) and <9 2 2 w„(fc). This is because for k along these rays we have (—ki,k 2 ) — r a (k) or (ki, —k 2 ) = r a (k), 
so that both (|2T5]) and (I2T9]) or (f2T20T) apply. In detail, suppose 

(-jfel, k 2 ) = r a (k), i.e. k = \k\e l{n/2 - a/2) or k = \k\ e -' l{ - n / 2+a /^ = ^kle^^- 01 ^ . 

Then it follows that 

d 2 kl uj n (k) = (dl 1 u) n )(-ki,k 2 ) = cos 2 (a)d 2 kl uj n (k) - sin{2a)d 2 klM uj n {k) + sin 2 (a)^ 2 w„(fc), 
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where the first equality holds due to (|2.19[) and the second due to (|2.15l) . As a result, for a G (— tt, tt], \a\ ^ tt/2, 
and k = ±\k\e^/ 2 - a / 2 ^ we get 

dl lM ^n{k) = itan(a) (d 2 u n (k) - S^u n {k)) . (2.22) 

Identity (|2.22l) applies also in the case when the S 2 reflection symmetry and the rotational symmetry (|2.12[) 
are both present for some a e (— it, — n], \a\ ^vr/2. Then (|2.22p holds for k that satisfy 

(ki,-k 2 ) = r a (k), i.e. k = ±\k\e~ ia/2 . 

2.3. Example: Hexagonal Lattice with a Circular Material Structure. As an example we con- 
sider the hexagonal lattice in the {x\, ir^-plane generated by the vectors 

a (1) = ao (s°n(I/3)) and fl(2) = a ° (0) witn °o > °- 

In the Wigner-Seitz cell U the material structure is given by the annulus centered at the lattice point in the 
origin and having outer and inner radii a$/2 and ao(1.31/4.9) respectively. The material properties are given 
by f](x) = 2.1025 for ao(l. 31/4.9) < \x\ < ao/2 and rj(x) = 1 otherwise. This is the same as the crystal 
used in [4] , where the corresponding band structure was also computed. One choice of vectors generating the 
reciprocal lattice is 

7,(1) — 22L ( a 2 2> \ _ 2tt / N 1.(2) _ 2w ( -a^ \ _ 2tt / 1 \ 
" ~ J12 I -a< 2) J ~ ao sin(ir/3) V 1 / ' U ~ J12 \ J ~ a \ ~ cot(7r/3) ) ) 

where J\ 2 = det(a^\ a' 2 ^) = a^af, 2 '' — a 2 af^. These vectors have been obtained via the formulas = 
27r s(^xiC3)_ and ^2) = 2 7r °i!^l« ) w herea« - (a« T ,0) T , = (6« T , 0) T for j e {1,2} and = (0,0, 1) T , 
cf. [5J Ch. 5] . Figure 12.11 shows the crystal geometry and the corresponding Brillouin zone. 

In this case both the rotational symmetry (12.12)) with a — tv/3, the reflection symmetry (I2.16|) as well 
as the analogous symmetry S 2 do hold. The band structure and Bloch waves can therefore be recovered 
via (j2~13l) . (j2~T4| and (pTT7| . (|2~T8)) from the irreducible Brillouin zone B in Figure [2TTT1 i.e. the triangle 
with vertices T,M,K, where T = (0,0) T , M = \b^ 2 \ and K = ^|& (2) | (1, 0) T . These points are called high 
symmetry points. 

Next we provide some specific information about the values of the second derivatives of u) n at the high 
symmetry points of the Brillouin zone at hand using symmetries (I2.15[) . (I2.19j) . and (|2.20p . This information 
will be used in Section I3T41 

Identity (|2.15p with k = and a = tt/3 yields 

d% 2 GJ„(T) = d 2 kl uj n {T) and d 2 klM u n {T) = for all n 6 N. (2.23) 
Symmetry (|2.20[) implies 

9l lM ^n(K) = for all n£N. (2.24) 
At k = r 2lr/3 (M) (= i& (1) ) we have k\ = so that ([2~T5]) implies 

9L^n(r 2l / 3 (M)) =0 for all n GN. (2.25) 

Relation (|2.22p then yields 



d 2 k2 oj n {r 27T/3 (Mj) = dlcj n (r 2w/3 (M)) for all n e N. 



(2.26) 
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(a) 

(b) 




Fig. 2.1: (a) Hexagonal lattice with a cylindrical material structure, (b) the corresponding first 
Brillouin zone B with a shaded irreducible Brillouin zone Bo. Note that the Brillouin zone has 
been scaled to fit the figure. 



Applying now (I2.15[) with a = 2tt/3, we get 

d 2 kl oj n (M) = d 2 ^ n (r 27l/3 (M)), d 2 k2 uj n (M) = G\>„(r 2V 3(M)), and d 2 klM u; n (M) = (2.27) 
for all n£N. Because r n / 3 (M) is obtained from r 27r /3(M) by the reflection (ki, fca) ~> {hi, — fea), we also have 
d 2 kl Lo n {M) = dlcj n (r w/3 (M)), d 2 k2 oj n (M) = d^ 1 u n (r T / 3 (M)), and d 2 lM uj n ( r7r/3 (M)) = 0. 

As an example we took the configuration from [3] as described in Section 12.31 The computations were 
done with a finite element Maxwell solver that uses lowest order Nedelec elements [22]. These elements were 
implemented in the software deal. II [7]. The eigenvalue problems were solved by a Krylov-Schur methodQ 

We computed the eigenvalues {k>n(&)}n=i,i4 and corresponding eigenfunctions {p n (k, •)}n=i,i4 f° r each 
vertex k in a discretization of the Brillouin zone B. The error level of this computations is about 10 -3 in the 
curl-norm and it is estimated from a series of computations on a sequence of nested grids. 

In Figure [2?2l we present the numerically computed band structure over OMq (following the tradition) for 
the above described crystal and for k = 5(27r/ao). Here, 9Bo is represented by 128 /c-points. It has, however, 
been checked that the observed gaps do not get narrower in the interior of B. Three band gaps appear on the 
positive half of the uj axis, one between and u)\, another one between uiq and 107 and the last one between 
W12 and W13. 

To get the extremal points at the band edges we used a bisection method in k which was initialised with 
data obtained from the band structure computation. The approximations to 1st and 2nd order derivatives of 
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k M> uj n (k) at the extremal values were obtained by first projecting k H> uj n (k) onto a locally quadratic finite 
element space and then taking mean values of the derivatives around vertices. 




Gamma M K Gamma 



Fig. 2.2: Band structure k i-> u n {k) for the described hexagonal lattice with the cylindrical 
material structure: the first 14 eigenvalues along 9Bo. Three band gaps appear on the positive 
half axis Q. > 0: one between and Wi, one between ojq and tJ7, and one between W12 and W13. 
Gap edges are marked by si, . . . , sg. 



3. Derivation of Coupled Mode Equations for Gap Solitons near Band Edges. 

3.1. Slowly varying envelope approach. We seek gap solitons E of (jl.9p . Afterward, the full electric 
field can be recovered via (|1.4I) . 

In the following let us assume that 
(Al) the spectrum {ui n (k) : k G B, n 6 N} possesses a gap, 

(A2) one of the gap edges, denoted by , is attained at precisely N 6 N points fc' 1 ^ , . . . , feW € B by bands 
with indices rii, . ..,njv, respectively, where the fc-points and/or band indices are not necessarily all 
distinct, 

(A3) for each j £ {1, ... , AT} the band uj nj is twice continuously differentiable in k at k — k^\ 
(A4) dlco nj (k^ ), the Hessian of w nj . at k = k^\ is (positive or negative) definite for each j E {1, . . . , N}. 
The smoothness assumption (A3) is needed to justify our Taylor expansions of u) nj near k^\ Bands ui n 
are generally only Lipschitz continuous due to possible transversal intersections of bands and their numbering 
according to size |20j . Away from points of intersection or tangency bands can be shown to be actually analytic 
in k by standard perturbation theory |18j . The simplest situation when (A3) is satisfied is thus when each 
band ui nj is isolated near k^\ which is equivalent to nx = . . . = ripf due to our ordering of bands according to 
size of w n (fc) at each k. 

Note that since each band w nj has an extremum at k = k^\ we have dk^nJk^') = dk 2 0J nj {k^) = for 
j £ {1, . . . , N}. Assumption (A4) then guarantees that the leading order terms in the Taylor expansion of the 
band w rLj around k — few) are in fact quadratic. 

The asymptotic expansion for the electrical field E of gap solitons with lu in the gap and in the vicinity 
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of the edge is expected [121 113] to be of the following slowly varying envelope form 

N 

e A o (»Kj 5 x) + e 2 ^ 1] (x) + e 3 ^ (2) {x)+0( 



3=1 

uj = + tte 2 , y = ex, < e < 1, 



(3.1) 



where Aj : M 2 — > C is a fast decaying smooth function and where Q = ±1. The sign of f2 is determined by the 
condition that + e 2 il lies in the gap. 

Performing a multiple scales analysis in the physical variables (x, y) is impossible. The reason is that in 
order to solve the resulting equations at each order of the expansion, one has to ensure that inhomogeneous 
terms are orthogonal to the kernel of L — 0^77, i.e., to u rlj (k^; . ) for all j e {1, . . . , N}. This orthogonality 
needs to be checked on the common period of those u nj . If, however, one of the components of few') is irrational, 
the corresponding u n , is not even periodic and this approach fails similarly to |13) . We therefore perform the 
asymptotic analysis in Bloch variables where all functions are [/-periodic in x and orthogonality conditions 
are always posed over U . 

Let us define the Bloch transform T : E 1— > E and its inverse, cf. [5J Ch. 7], by 

E(k;x) = (TE)(k;x)= ^ e iK - x E(k + K), E(x) = (T^ 1 E)(x) = j e ifc ' x E(fe; x) dfc 

for all x, k £ R 2 , where E denotes the Fourier transform of E 

E(k) := (.F£)(fe) := 7-^- / E(x)e~ ik - X dx, E{x) = {T^E^k) := [ E(k)e ik - X dfc. 



(271") 2 j - .1 - 

By definition we have the following properties of the Bloch transform 

E(k; x + R) = E(k; x) for all Re A, 

E(k + K; x) = e~ lK x E(k] x) for all K £ A* . (3.2) 
Multiplication of two functions /, g in physical space corresponds to convolution in Bloch space, i.e., 



(T(fg))(k;x)= / f(k-l;x)g(l;x)dl=: (/ * m g)(k;x), 

where (|3.2p is used if k — I £ B. Especially, if x v- > f(x) is [/-periodic, then 

(T(fg))(k;x) = f(x)(Tg){k;x). 

This can be easily checked by writing / in the form of a Fourier series, i.e. f(x) = Y^KeA* c k^ iKx , cf. [8j 
Ch. 7]. Exploiting this observation and applying the Bloch transform to (|1.9|l leads to 

(Z(k) -uj 2 n(x)J E(k; x) — lj PNL(fc; x) for all x, k £ R , 

where 

^NL(fe; •) = X i I ) T{2\E\ 2 E + E-EE) = X ^(2(£.* B £)*»£ + (£.*,£) * m E), 

with / . * B g := J2j fj *m 9j f° r vector valued /, g, while / * B g is understood componentwise for scalar / and 
vector valued g. By definition of the Bloch- and Fourier transformation one immediately finds 

r(A J (e.)e ikU) <^)(k- 1 x) =e- 2 £ Aj ^(fc - fc u) + K)^j e iKx , 

KGA* 
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so that the asymptotic ansatz (|3.ip is transformed to 

s4> {0) (k;x) + e 2 4> (1) (k;x) + e 3 4> (2) (k;x) + 0(e 4 ), (3.3) 

where 

JV 

ft°)(k;x) =e- 2 J2 53 4* (H k - kU) + K )) e iK - x Pn 3 (k Uh ,x). 
j=i KeA* 

Similarly to [12] and [T3], due to the fast decay of the Bloch transform of Aj in k, we approximate Aj (j(k — few + 
by Xb,t (& — & + ^0 A? (|(& — ^ + ^0) f° r some r G (0, 1), where \S is the indicator function of a set S 1 , 
D 5 :=B S (0) with B 5 (z) := {fc e M 2 : |fe - z| < (5} for J > 0, z G M 2 . 
We will therefore introduce the approximation 

£(fe; a?) = e- 1 ^ ' (fc; x) + (fc; x) + (fe; x) + 0{e 2 ) 

with 

AT 

&°\k;x) = 53 53 *JV ( fc _ fe0) + K ) ^ (e( fc _ fc0) + *0) e^^Pn^fe 00 ;^ 
i=i KeA* 

for all fc G B and i£l 2 . In the following we will use the notation K m — mib^ + ni2b^ G A* for m G Z 2 for 
convenience. As an abbreviation we let l^' m \k) := |(fc - fcW + if" 1 ) for fc G R 2 and m G Z 2 , so that Z5 (0) is 
given as 

N 

3°\k;x) =53 E X DE 4^ m Hfc))^(^ m Hfc))e iKm -^ 3 (fc°^)- (3.4) 

Note that £ ,< -°- > ( . ; a;) is supported on a set of (for sufficiently small e) disjoint balls B e r{k^ — K m ), j G 
{1,...,N}, m G Z 2 . 

3.2. Formal asymptotic analysis. Let us proceed with a formal asymptotic analysis of (|1.9j) . First, 
we consider fc close to few" - K m , i.e., fc G B £ r (fe^ - K nl ) for some j G {1, . . . , N}, meZ 2 . Then 

Z(fc) = Z(fc^ - K m + £l^' ro )(fe)) 

= Z(fc w - K m ) + e^ m) (fc) ■ d k L{k {]) - K m ) + ie 2 Q(^ Vm) (fc)), 

where we have used the fact that the second derivatives of L are constant in fc, see (12.81) . and where 

2 

^• m) (fc) • d k L(k^ - K m ) =Y,4 j ' m) (k)d ki L(k ( - j) - K m ), and 

i=l 
2 

Q(^ m \k))= 53 ^ m) (fc)^ m) (fc)^ 2 a , fei) z. 

a, 6=1 

Using (I3.3[) , p. 41) , (|3.5p and to = cu* + He 2 , we get a hierarchy of equations at each power of e for x G U 
and fc G B E r{k { ri - K m ). We now study the equations related to e 1 ,e°,e 1 under the condition that the 
nonlinear term contributes to e , which is confirmed later in (I3.14[) . 

0(e _1 ): The resulting equation is 

A^l^ik)) (Z(fc« - K m ) - u£r}(xf) (p nj (k^;x)e iKm - x ) 

= Aj(£ U ' m) (k))e iKm - x (Z(fc w ) - c#j(a:)) Pnj {k (3) -x) = 0. 
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This holds by the definitions of lj* — Lo nj (k^) and Pnj (few ; . ). 
O(l): The resulting equation is 

(Z(kW - K m ) - uln{x)) E^(k;x) 

= -A 3 {t^ m \k))(t^ m \k) -d k L(k^ ~K m ))( Pnj (k^;x)e iKm - x ) 
= -Aj(^«- m )(fc))e ix " , - a; (^ m )(fc) • d k L{k^))p nj {k^;x) = 0. 

Using (|2.6p . the solution is found to be 

E^(k-x) =Aj(^(k))e iKm - x (£^ m \k) -d kPnj (k^;x)), (3.6) 

where (fc) . d kPn } (k « ; x) = lf m) (k)d kiPrij (k^;x). 
O(e): The contribution of L(k)E is 

(Z(fc( J ') - if m ) - a#?(s)) & 2 \k-x) 

+ ±Q(^' m) )£ (0) (fc; z) - 2w,nj7(as)^°5 (fc; x) 
+ (t^' m \k) ■ d k L{k {i) - K m ))E {l \k- lX ). 

By insertion of the previous results this gives (for k G B £ r(k^ - K m )) 

(Z{k® ~ K m ) - a#j(z)) E (2) (k;x) 

+ [^Q(^' m \k))p nj (k^;x)-2^n V (x)p nj (k^-,x) 

+ (£^- m) (k) ■ d k L(k U) - K m )) (£^- m) (k) ■ d kPn] (fc (j) ; a:))] A 3 (£^- m) {k))e lI<m - x 

= ^ 2 X^W^(2(^ (0) .*,fM) * B £ (0) + (£ (0) .*. £ (0) ) (3.7) 
=:ulxf{x)G 3 ik,x). 

The remainder of the section is devoted to the analysis of the structure of Gj in (|3 . 7[) and to the derivation 
of a solvability condition for Q3.7P . 

Let us first analyze the nonlinearity. The convolutions in (|3.7[) can be expanded into the form 



A' 



a,/3,7=l 

where a, 6, c e {1, 2, 3}, and functions £ Q . a and £*_ a are given by 

Ua{k; x) := p„ QiQ (fc( Q ); x) £ XnA k ~ fc(a) + (l( fc ~ fcW + eUS *'*> 

*) :=|^(fc( Q );x) £ x DEr (fc + fc(«)-^)2 Q (i(fc + fc( Q )-^))e- i ^-. 

zez 2 

Note that (13.81) represents all the nonlinear terms in (|3.7[) due to commutativity of *„. The summands in 
have the form 



(3.8) 



(?a,a*i^,l*B^ )C )(fc;a;)= X! 9noq{k)X)= ^ 



9noq\k] x), 



(3.9) 



n.o,q^l 
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where (with indices a, /3, 7, a, b, c suppressed) 

9noq{k',x) ^ Pn a ,a{k^ ^ x)p n/3 {kS^ ] x)p n _ J t c(k^ x) 

X Ds r (k-t- fc (a) + K n )A a - 1 - fcW + K n )) 

x X De r (t-s- + X°)A,3 - s - jfeW + K°) 
x XD er (s + — K q yA 1 (i(s + fc (7) - if 9 )) dsdi 

and with 

M 7 = {z e Z 2 : fc - k {l) + K z E B e r(0) for some k E B and all e > 0}, 
M, (2) = {z G Z 2 : fc - fc (b) +K Z £ B e r(0) for some fc £ B + B and all e > 0} 



(3.10) 



for b G { a , f3 } . The truncation of the series in (13.9[) comes from the fact that for s , t , k £ B we have t—s£ 

v 



and fc— t £ B+B so that the three characteristic functions in (|3.10l) can be nonzero only for n £ M a , o £ M\ , 



and q £ My. More precisely, this is seen as follows. 

Only those combinations of n, o, q which produce nonzero values of all the three characteristic functions in 
(|3.10P and of the function Xo e r ( ■ — & +K m ) in (|3.T[) for given j and some k,t,s EM are of relevance. Firstly, 
Xd t ( s + ^ c< ' 7 ' ) — K q ) i s nonzero for some s £ B and for arbitrary e > if and only if so := — k^ + K q £ B 
(the closure of B) for some q El?, which is equivalent to 

q E My. (3.11) 

Secondly, for a fixed q the factor %£, r (t — s — k^> + K°) is nonzero for all e > and some f £l and s obtained 
in the first step if and only if to '■= sq + k^' — K° £ B, i.e., 

k (f3) - fc (7) + K q - K E B. (3.12) 

This can always be satisfied by a choice of o £ M^ . Finally, for fixed q and o we need that \d r{k — t — k^ + 
K n ) does not vanish for some k £ B with k — fcW + K m E D e r and all e > 0, where this latter restriction is 
due to the restriction k £ B s r{k^' — K m ) in (|3 . 7[) . In other words, we need that ko := k^' — K m £ B and 
= fc - to - k^ +K n , i.e. 

fcCo) + _ fcW + _ _ if" = fcW) - K m g I (3.13) 

for some neZ 2 . In fact, all solutions for n of (|3.13|) lie in M a 2 \ 

In summary, for a, /3, 7 G {1, . . . , iV} the term g noq is nonzero in (|3.10p if n, o, q satisfy (|3.1ip . (I3.12[) and 
(|3.13[) . So the term £ aiQ * B * B enters Gj provided 

A a ,p, 7}j := |(n,o,g) G (Z 2 ) 3 : n G M a 2 \ o £ Mf \ q £ M 7 and (l3~T2l) . (l3T3ll holdj 

is nonempty for some m E Mj. Note that we omitted an index m in this definition, because the set is either 
nonempty or empty for all m £ Mj. Indeed, if m is one index that meets the requirements with (n, o, q) and 
z is any other index in Mj, then z meets the requirements for (n + m — z, o, q). Aa j can be constructed 
by a computer code that scans all possible combinations of n,o,q. This will be discussed in Section [3.41 

Due to the characteristic function, the integration domains in f|3 . 10[) can be reduced to s £ B e r(—k^' + 
K q )r\M&ndt E B 2s r(k^-k^-K +K q )nM. Now we introduce the change of variables s := (s+k^-K q )/e 
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and t := {t - feW + JfeW + K° - K") /e to get 

9noq{k\%) —E e ^ Pn a ,a(,k ^ ! x )Pnfj .b{k^^ x)Pn~, ,c{k ^ j 2<) 

Xd r 



i>„ r „in Ur-in 



fc-(fc ( ° ) +fc ( ' 3) -fc (7) )+i<'"+i<'°-j<" i _ 

(3.14) 



x £ ( fc - (fc( ° )+fc( "- fc( ; ))+ "" + "°-" g - t) Xn^ (i ~ s)Mi ~ ~s)Xn er _ 1 ®% (S) dS dt. 

The factor e 4 in this formula shows that Gj — 0(1) as required for the consistent asymptotic expansion. If 
(|3.13l) is satisfied, (13.141) becomes 

j j x D ^(^^-t) 



(3.15) 

B+fctT) -K1 



a ( K k - ki3) + Km - t) Xd^ {t-s)Mt-~s)x Der _ 1 (5)^(5) dSdf 



(3.16) 



x A 

for fc G B e r{k^> — K m ). As we show in Remark 13.11 summing, for fixed k,j,m, the terms (|3. 151) over 
(n,o,q) € A ai fj n j yields a double convolution integral in s,t over the full discs D 2e r-i and D e r~i, i.e., 

*£, c )(fc;s)=eV<* Co, +* W - fcW -^ 

= . e 4 e i(- fe O) + ^)., Un ^ a(fc ( Q ). x)u ^ fc(fc (/ 3);x)lI _ (A .( 7 ). a;) ^^(^(fe)) 

for £ B e r(k^ s> ~ K m ). Here we have used u na ^ a (k^; x) — p na ,a(k^ ', x)e ika ' x , etc., and we defined 
K ryi^'^ik)) as an abbreviation for the integral on the right hand side. 

Remark 3.1. To show that the sum of g noq over (n,o,q) G A a ,/3,-yj yields a double convolution integral 

(2) 

over full discs, let us first note that the definitions of M 1 and ensure 

(J ((B + fcW - K*) n iv) = IV, (3.17) 



|J ((B - kW + fc (7) + #° - K q ) n £ 2 e") = ^2 £ " • (3.18) 

These are obvious when feW G int(B) and k^\ fcM € int(B), respectively, because thenM 7 = Mf ] = {(0,0) T }. 
But when k^' G 9B, i/ien on/y a fraction of —k^ + D £ r lies in B fin our example with a hexagonal B the 
fraction is a half unless k^ is a vertex ofM, in which case it is a third) and the rest lies in periodicity cells 
centered at neighboring reciprocal lattice points. Each point i in this rest is therefore mapped to B via £+K q with 
some q € Mj, and we thus have (|3 . 1 7|) . By an analogous argument, observing that k^ — (k^ — K q ) G B + D e r 
for all q G M 7 , we get (|3.18|) from the definition of Mo . 

Let us now assume (|3.13[) and show that for each K q fixed, i.e. for each fixed integration domain in the 
inner integral in (|3.14p . the sum of g noq over (n,o,q) G A a ,p,-y,j yields an integration over the full disc D 2e ^ 
in the outer integral. If this were not the case, i.e. if 31 G D 2e r such that £ <£ B - fcW + fcM + K° - K q 
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for any such (n,o,q) G A a jj.yj, then by (|3.18[) there would be o G Mo such that (n,o, g) A a ,p^,j while 
(|3.12p and (|3.13|) are satisfied. This is a contradiction to the definition of A a ,p^,j ■ After that we sum over 
all q G My and the result follows from (|3.14l) . 

We now write the d-th component (d G {1, 2, 3}) of Gj as 

G 3 Ah, x) = e-^XDer (k fcCfl + K m ) J2 r ilc fe 0) ^ 0) ^F) (*5 (3-19) 

aM,c=l \ ' 

where the integer coefficients T^l can be easily derived from (|1.7|) . In detail we have 1 = T^l 2 = Tg 3 ^ 3 = 

o r (l) _ r (l) _ r (l) _ P (lj ' _ r (l) _i r (2) _ r (2) _ r (2) _p(2) 'l r (2) ' '_ r (2) '' 

°> 1 1,2,2 — 1 2,1,2 — 1 1,3,3 — 1 3,1,3 ~~ 1 2,2,1 — 1 3,3,1 ~~ x ' 1 1,2,1 — 1 2,1,1 — 1 3,2,3 — 1 2,3,3 — 1 1,1,2 — 1 3,3,2 — x i 

rS,i = rg,i = rg j2 = rg >2 = r^ 3 = rg )3 = 1, and the remaining rJ$ )C are zero. Finally, using (031, 
we get for jfe G B s r(k^ - K m ) 

o,6,c=l a,/3, 7 e{l,...,JV} s.t. (3.20) 

xu^(fc (7) ;a;)ftS, 7 (^'' m) (*))- 

In order to make the discussion of the asymptotic hierarchy complete, we also have to consider the part of 
the A;— domain outside the neig hborhoods of k {] \ For k G B such that k - jfeW + IT 71 el\D £r for all m G Mj 
we have - iT m ; 2) - w 2 rj(a;))EW(fc; a;) = for / G {0,1} so that E^(k; .) = E^(k; .) = for such k. 



3.3. Coupled mode equations. We return now to equation (|3.7p . Due to the Fredholm alterna- 
tive the existence of A-periodic solutions of equation (13. 7[) is equivalent to L 2 -orthogonality of (|3.7|) to 
p nj (k^;x)e' K ' x , which needs to be ensured for all m G Mj and j G {1, . . . ,N}. The range of ^C?' m ) is a 
different section of the disc D e r-i for each m G Mj. This section is a (l/|Mj|)-th of the full disc so that these 
\Mj\ equations actually build one equation in I G D e r-i. Figure [3~T1 shows these sections for two example 
points few. One example is for \Mj\ = 2 and the other one for \Mj \ = 3. 

When imposing the orthogonality condition, the common factor e lRm ' x of the right hand side of (|3.7[) is 
canceled in the complex inner product with p n . (k^ ; x)e lK ' x , so that the same solvability condition holds for 
all m G Mj. Using (|3.7[) . (|2.5p . and (|2.7p (with and fc + replaced by rij and fc^), we obtain 



QAjy 



- \ + ^ 2 w„, (*<*>) + n^d\ uk ^ ni 2, (£) + A/} W = (3.21) 

for I G D e r-i , where 

^•W = Y<X^(-)Gi(£; .),fMfc«; .) e ^ m '(-)) 



T E r Sc E / x^W^, Q (fc (Q) ;x) Uri „ fc (fc(«;x) 

x ti^(feW;a;)u^5(A;W;a;)dc ^ i7 (^) 



,&,c,rf=l a,/3,7G{l,...,iV} s.t. 



a,/3,7e{l,...,JV} s.t. 
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(a) 



(b) 




Ber(*W)) 




m^iK fed) 

B e r{k^) 



Fig. 3.1: Two example points k^' in the case of the hexagonal lattice and the corresponding 
ranges of e£ u ' m) for all m £ Mj. In (a) we have Mj = {(0, 0) T , (1, 1) T } =: {m (1) ,m (2) } and 
in (b) Mj = {(0,0) T , (0, 1) T , (1, 1) T } =: {m (1) , m (2) , m (3) }. The shaded sections along the 
boundary of B are those k £ B for which xc e r (k — + K m ) 7^ for the m £ Mj written next 
to the respective section. 



i.e. with (fl~7) and the definition of V in (|3~19| 



a,b,c,d— 1 



T 



[2K a (fcW ; -)K,(fc ( «; .) 

+KJfc( Q ) ; .)- u „,(^) ; .))«^"(* (7) ;-)] -^7(fc (j) ; ■) 



(3.22) 



The symmetries in rj, d ? imply symmetries in l a ,0,<y,j- Namely, due to the symmetries It9 = and 
rile = rS, rf we have' ' 



^a,/3,7,j =Ip,an,i and ^o,/3,7,J = -W.j/y for all a, /3, 7, j £ {1, ...,iV}, 



-.(d) _ p(6) 
a.b.c c.d.a 



Ia,p,i,j = Ij,j.a,p for all a,/3,7,j £ {1, . . . ,iV}. 



(3.23) 



(3.24) 



Symmetries p.23[) and (|3.24[) imply, in particular, that I a ,p, a ,p — I a ,p,p,a G K f° r au a, /? £ {1, • • • , AT}. 

Let the crystal satisfy the rotational symmetry 77(2;) = -q{r v {x)) and xii O 2 ') = Xci^( r i'( a; )) f° r ai l 33 G ^ 2 
and some v £ (— 7r, 7t] and let U be chosen so that r v (U) = U. If for each m £ {a, (3, 7, j} C {1, . . . , iV} there 
exists m! £ {1, . . . , N} such that 

fc( m ') = 7v(/fc( m )), 

and if ^(/c^™- 1 ) is a geometrically simple eigenvalue of (|2.3p for each m £ {a,/3,7,j}, then 

^a,/3,7,j — la' ,/3 f ,*y',j f • (3.25) 
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This is seen by the change of variables y = r u (x) in (|3.22p . using the facts r v {U) = U and r v [v) ■ r v {w) — v ■ w 
for all v,w £ C 3 , and employing the symmetry (I2.14[) . 

(3) 

Additional symmetries in I a ,p,jj arise when a spatial reflection symmetry in r\ and xli i s present. For 
instance when rj(x) = n(S2{x)), Xci 0*0 = X^\^2{x)) for all x £ M 2 (see Section |2.2.2[) and if for each 
m e {a, f3,j,j} C {1, . . . , N} there exists m' £ {1, . . . , N} such that 

fcK) = S 2 (k {m ' ) ) 

and such that S , 2(^ m - ) ) = fc( m ) does not hold for any m £ {a, /3,7, j}, then 

Ia,P,"/,j = Ict',0',7',j'- (3.26) 

This is proved via a change of variables in (I3.22|) using (|2.2ip . where a = due to our assumptions. A similar 
result holds for the reflection symmetry in X\. 

Returning now back to (|3.2ip . for smooth envelopes Aj we can neglect the contribution of Aj from 
I € M 2 \ D E r-i or simply assume that Aj satisfy (13.21)) also there. This step can be rigorously justified via a 

persistence argument similar to that in [13l [14]. h^p 7 will then be replaced by A a * Ap * A 7 . The inverse 
Fourier transform then produces the coupled mode equations 

QAj + \ (dlu n] (k^)d 2 yi + d 2 k2 u nj (k^)d 2 y2 + 2d 2 kiM u nj (k^)d 2 yuy2 ) A, +Afj = (3.27) 

on R 2 , where Afj is given by 

,/3,J ,jA a ApAry. 

a,/3,~f£{l,...,N} s.t. 

Note that the coupled mode equations have the same general structure as those for gap solitons of the scalar 
Gross-Pitaevskii equation [TJ] . 

A localized solution A of (|3.27[) should produce via (|3.1j) an approximation of a gap soliton of the Maxwell 
problem (|1.9[) . A rigorous justification of this statement can be done via the Lyapunov-Schmidt reduction 
similarly to [121 E3 HI] and will be the subject of a future project. System (13.27)) does not have localized 
solutions for arbitrary values of coefficients. The coefficients of the derivative terms are given by the band 
structure and CI = ±1 is determined by the condition that w = lj* + s 2 Vt lies in the gap. But the function yffl 
in I a p ~ j has not been fixed and remains free at this point. 

The linear part of the operator in (|3.27p is definite due to our assumption (A4) in Section [3TT1 and the fact 
that Cl < at upper edges and f2 > at lower edges. The linear part of the operator is positive definite at 
lower edges w*, where k^' are points of maxima and negative definite at upper edges. In case N = 1, where 
Ml = 7|A 1 | 2 A 1 and 7 = SuX^i\ u nj{k^\ -)| 4 , a localized solution exists in the upper edge case only if 

(3) (3) 

X c i is such that 7 > while in the lower edge case x c i has to produce 7 < 0. Physically it makes sense to set 

Xci = there, where 77 = 1 (i.e. in vacuum/air). In the annulus regions, where i] = 2.1025, we set x C i = 1 

(a focusing nonlinearity) if 7 > is needed and x ci = — 1 (a dcfocusing nonlinearity) if 7 < is required. 

This is in agreement with previous results on bifurcation of gap solitons from spectral edges in the periodic 

nonlinear Schrodinger equation [T71 131 H31 HH where bifurcation from upper/lower edges occurs for the 

focusing/defocusing nonlinearity respectively. In the case N > 1 our numerical examples produce all I a ,p.~t,j 

(3) 

of the same sign so that we set in the annulus regions, once again, x c \ = 1 if is an upper edge of a gap 
and Xti = —1 if it is a lower edge. 



OUT-OF-PLANE GAP SOLITONS IN 2D PHOTONIC CRYSTALS 



19 



3.4. Examples of Coupled Mode Equations. We present next coupled mode equations for gap 
solitons in the vicinity of spectral edges for the example in Section [2. 31 as well as for other canonical examples. 
As seen in Figure [2T2l there are 3 spectral gaps (0, si), (s2, S3) and (S4, S5) on the positive part of the spectral 
u> axis for this specific example. We have the numerical values 

si = wi(r) fa 3.610, s 2 = w 6 (r) fa 3.701, s 3 = w 7 (r) ps 3.750, 
s 4 = wi 2 (0, 2.351) ps 3.873, s 5 = wi 3 (0, 2.407) « 3.882. 

At sx and S2 several bands lie very close to each other at the extremal point k = T. It is, however, not known 
whether these truly touch and the eigenvalues have higher multiplicity than one. Numerical tests have shown 
that varying the value of rj for the annulus material does not change the ordering of bands at k = T near Si 
and S2- We, therefore, assume that the edges si and S2 are simple eigenvalues at k = T leading to N = 1 at 
si and s 2 . If it can be proved that, for instance, s\ is indeed a double eigenvalue, then N = 2 at s\. Likewise, 
N would change if the multiplicity could be established for S2 ■ 

Similarly, the band uj\2 is close to u> = S5 at four distinct fc-points along 8Mq. At the point k = (0, 2.351) 
the numerical value is maximal and an analogous test shows that it remains maximal for a range of values of 
77. We thus assume that within Bo the value w = S5 is attained only at k = (0,2.351). Due to the discrete 
rotational symmetry of the band structure we thus have N = 6 at S5. Analogously, we have N = 6 at S4. 

Except for the simplest case with N = 1, like in Section T3.4.11 we determine the sets A a ,p n ,j using 
a Matlab program. First of all, it is clear that for any fc b £ B the sets M\, and contain only those 

(n,o,q) £ Z 2 with m,oi,qi £ { — 1,0,1} for I — 1,2. To determine A a ,/3,-y,j, we therefore need to test only 
finitely many integer vectors (n,o,q) for conditions p,12[) . (I3.13[) . For an example with N = 3 we show in 
Section [3.4.21 the resulting sets A a .p^.j computed using this routine. 

3.4.1. Coupled Mode Equations near Edges for the Example in Section 12.31 

Coupled Mode Equations near the Edges s\, S2 and S3 (N = 1). 

At the edges s\, S2 and S3 in Figure [2T2l the situation is particularly simple. As discussed at the beginning 
of Section EH we have N = 1 and fcW = T = (g). Since £ int(B), any (small) neighborhood of k^ 
lies completely within B and thus Mi = {(q)}- A simple inspection determines that we have -4i, 1.1,1 = 
{((o)>(o)>(8))}' T ne resu h m g coupled mode equation for A = A\ is 

(n + a(dl+d 2 y2 ))A + 1 \A\ 2 A = 0, (3.28) 

where a = \d 2 ki uj ni {T) = \d 2 k2 Lu ni {T) (cf. O) and 7 = h, lxl . 

The three cases si, S2 and s 3 differ by the value of ni, i.e. the band index. At = s± we have n\ = 1, 
at = S2 we have n\ = 6 and at = s 3 we have n\ — 7. And, as discussed at the end of Section 13.31 at 
the upper edges si, s 3 we have fi = — 1 and the function %c? nas the value 1 in the annulus regions and 

f 31 

otherwise. At S2 we have f2 = 1 and x c { = — 1 in the annuli. 

In Section 14.11 we present a numerical example on a gap soliton approximation near S2. We list here, 
therefore, the numerical values of the CME coefficients for the case S2 : 

w« = s 2 ~ 3.701 : a « -0.0107, 7 w -3.057. 

Coupled Mode Equations near the Edge S5 (N = 6). 

At the upper edge S5 in Figure [2~2l we have N — 6, n\ — n2 — ■ ■ ■ = rig = 13, k^ pa (0, 2.458) lying on 
the line from T to r 2 ^/ 3 (M), and k^\ j = 2, . . . , 6, obtained via a rotation of fc^. In detail 

k U) =r, , 1)2 r(fc (1) ) fori = 2,..., 6. 
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The symmetry properties (I2.15|) . (|2.19p . and (|2.20[) produce relations among the linear coefficients of the 

CMEs. The sets Aa^^j are either empty or contain solely the element {(o))(o)>(o)} as checked by the 
Matlab routine. The resulting CMEs are 

(n + ai8^+p 1 ^ a )A 1 + N 1 = 0, 

(n + a 2 dl + p 2 d 2 y2 + nSfa m )A 2 + N 2 = 0, 

(n + a 2 dl+f3 2 dl-vdl, y2 )A 3 + N 3 = 0, 

(n + ai^ i + ( 8 1 ^ 2 )A 4 + iV4= 0, 

(SI + a 2 d 2 yi + fi 2 8 2 y2 + fid 2 yi m )A 5 + N 5 = 0, 

(SI + a 2 d 2 yi + /3 2 d 2 2 - fid yi , y2 )A 6 + N 6 = 0, 

where fi = -1, <*i = ^^(feW), A = ^^(fcW), a 2 = i(a x + 3/3i), ft = ±(3c*i +&), M = ^K-ft) = 
C fc2 ^i3(fc (2) ), and 

6 

A^i = 253/ <iMi i|Ai| a Ai -/i,i,i,i|Ai| 2 Ai + 2(/ 2:5 ,4,i^2^5 +/ 3 ,m^3A)A 4 , 

i=l 
6 

N 2 = 2^2l t ^ l . 2 \A l \ 2 A 2 - I 2a , 2 , 2 \A 2 \ 2 A 2 + 2(71,4,5,2^1^4 + 73,6,5,2^3 A>)^5, 
»=i 

6 

N 3 = 2 ^ J 4A4 ,s|A | 2 A 3 - 7 3 , 3 , 3 , 3 |A 3 | 2 yl 3 + 2(7i, 4 ,6,3AlA 4 + 72,5,6,3^5)^6, 
i=l 
6 

7V 4 = 2^ 7 4 ,4, 4 , 4 |A 4 | 2 A 4 -7 4 ,4,4,4|^4| 2 A 4 +2(/ 2 ,5,l,4^2^5 + 4,6,1,4^M 6 )^1 , 
»=1 
6 

N 5 = 2j2k5^,5\A l \ 2 A 5 - 75,5,5,5 1 A 5 | 2 yl5 + 2(7i, 4 , 2 , 5 Ai A 4 + 7 3 , 6 ,2, 5 ^3^6)^2, 
i=i 

6 

6 - 7 6 , 6 ,6, 6 |A 6 | 2 A 6 + 2(7i, 4 , 3 , 6 ^i^l 4 + 72,5,3,6^2^5)^3. 

i=l 

Due to symmetries, many of the coefficients in the nonlinear terms are equal. Symmetry p.25[) with v = tt/3 
and symmetry (|3 . 23[) imply 



Using (|3~23|) and (|3~2"3l) . 



7o 


:=7i, 1,1,1 


= 4,2,2,2 


= • • • = 76,6,6,6, 




7i 


: =72, 1,2,1 


= 4,2,3,2 


= • ■ • = 7e,5,6,5 = 4,6,1,6 






=7l,2,l,2 


= 4,3,2,3 


= • • • = 7s, 6, 5, 6 = 76,1,6,1, 




72 


: =4, 1,3,1 


= 4,2,4,2 


= 75,3,5,3 = 4,4,6,4 = 4,5,1,5 


= 4,6,2,6 




=4,3,1,3 


= 72,4,2,4 


= 7 3 ,5, 3 ,5 = 74,6,4,6 = 4,1,5,1 


= 4,2,6,2 


73 


: ~ 4,1,4,1 


= 4,2,5,2 


= 76, 3 ,6,3 = 4,4,1,4 = 72,5,2,5 


= 7 3 ,6,3,6 


74 


:= 4,5,1,4 


= 7 3 ,6,2,5 


= 4,4,3,6 = 4,5,4,1 = 7 3 ,6,5,2 


= 4,4,6,3 


we 


get 








74 


= 73,6,4,1 


= 73,6,1,4 


= 4,4,5,2 = 4,4,2,5 = 7 2 ,5,6,3 


= 72,5,3,6 



We have 70 , 7i , 72 , 73 G K as explained below (|3.24l) . 
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Finally, because feW = (k[ 4) , -k ( 2 4) ) T , k^ = (fcf \ -fc< 3) ) T , and k^ = {k{ 6) ,-kf ] ) T with fc«,fc( 2 ) and 
fc( 5 ) lying in the interior of B away from the line k2 = 0, the symmetry p.26[) applies and we get 

^2,5,4,1 = -^3,6,1,4- 

Therefore 



12, 5,4,1 — ^3,6,1,4 — -^3,6,4,1 — I4, 1,5, 2 — ^2,5,4,1 (3.30) 

so that also 74 e M. The second, third and fourth equalities in (|3.30|l hold due to (|3.23|) . ()3.25j) . and (|3.24j) . 
As a result the nonlinear terms in (|3.29[) can be simplified to 



At 


:=2 


;>r 


|-7i(^2| 2 - 


HAI 2 )- 


f-72(^ 3 | 2 - 


HAI 2 )- 


1- T3 1 -4_4 1 2 ^ 


Ai- 


h 2j 4 (A 2 A 5 - 


1- ^ 3 A 6 )^I 


N 2 


:=2 




l-7i(^i| 2 - 


H^l 2 )- 


f- T2 ( 1 2 - 


HAI 2 )- 




A 2 - 


h 274(^1^4 - 


hA 3 Ao]A 5 


A 3 


:=2 




l-7i(^2| 2 - 


M^l 2 )- 


H72(-4i| 2 - 


HAI 2 )- 




A 3 - 


h 27 4 (AiAi - 


V A 2 A 5 )A^ 


A 4 


:=2 




l-7i(^3| 2 - 


HA| 2 )- 


l"72( ^2| 2 - 


H^l 2 )- 


^73|-4i| 2 J 


Ai- 


h 274(^2^5 - 


\- A 3 A 6 )A^ 


A 5 


:=2 




H7i(^ 4 | 2 - 


HAI 2 )- 


l-72(v4i| 2 - 


HV)- 


^73^2^ 




h 274(^4^4 - 


VA 3 A 6 )~A 2 


A 6 


:=2 




f 7i(^i| 2 " 


HAI 2 )- 


f- 72( -4 2 | 2 - 


H^l 2 )- 




A 6 - 


h 274(^4^4 - 


1- A 2 A 5 )A^ 



with 70, 71, 72, 73, 74 6 R. A system of six CMEs with the same structure as above arises also at the edge S4. 

In Section 14.21 a numerical example of gap soliton asymptotics near S5 is given. The numerical values of 
the coefficients in the CMEs (|3.29p for S5 are 

= s 5 ~ 3.882 : a x w 0.0189, a 2 « 0.146, ft « 0.189, /3 2 w 0.0614, ^ w -0.0736, 
70 « 1.282, 71 « 0.789,7 2 ~ 0.757, 73 ~ 1.193, 74 « 0.714. 

( 3^ 

As S5 is an upper edge edge, the coefficients 7?, J G {0, ... ,4}, were computed using = 1 in the annulus 
regions. 

3.4.2. Additional CME Examples. 

Example of Coupled Mode Equations for N = 2. 

An example of a situation for N — 2 is when the locations of the extrema are k^ = K, k^ = r^/ 3 (K). 

With &W,fe( 2 ) as in Section l2~3l we then have = ^ ( J) and k^ = ^ The corresponding integer 

shift sets are Mi = {(q ) i ( 1 ) i ( 1 )}) ^2 = {(0 ) > ( ) > ( 1 )}• Due *° * ne rotation symmetry of the bands and 
their labeling according to size, we necessarily have n\ = n 2 . We define n* := ri\ = n 2 . From <|2.24|) we have 

d 2 kl , k2 uj n ,(k^) = 

and using (|2.15p with a = 7r/3, wc obtain 

dlu n ,{k {2) ) = H^».(* (1) )+3^«>».(* (1) )), 

dl^Ak^) = i(3flg l «».(*W) + flg,^.(*W)) I 

^, fc2 ^,(fc (2) ) = f (^^(feW) -^(JbW)). 

After having numerically checked the sets A a ,p,y,j for all combinations of a., ft, 7,j to determine the nonlinear 
terms, we thus arrive at the CMEs 

(ft + ai c£ + ft^jAx + (70 Ail 2 + 2 7l |^ 2 | 2 )^i = 0, 
(ft + a 2 c£ + (3 2 d 2 y2 + »d* uy2 )A 2 + ( l0 \A 2 \ 2 + 2 7l |A 1 | 2 )A 2 = 0, 
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ft 



where ttl = ^^.(fcW), ft = ^w^kW), and a 2 = i(ai + 3ft), ft = ±(3ai + ft), fi 
7o := ^1,1,1,1 = h,2,2a using symmetry (|3.25p with v — n/3, and 71 := 2i,2,i,2 = ^1,2,2,1 using (|3.23|l . 
Example of Coupled Mode Equations for N = 3.. 

Let us assume that a gap edge for N = 3 has extremal points at — M, — r n / 3 (M), k^ = 
r 2 ^/ 3 (M). With the choice of the reciprocal lattice vectors V- l \b^> as in Section [212 we have fc« = ±M 2 \ 
fc( 2 ) = ±(&W +fe (2) ), and = with the corresponding integer shift sets Mi = {(§),(?)}, M 2 = 

{(0) >(!)}) an< ^ -^3 = {(0) ' (o)}- Similarly to Section [3.4.21 we have ri\ = n% = n 3 = : n*. Using (|2.15[) and 
(p ^l - fOTjl . we get 



The sets ^ a)/ g l7) j are, once again, determined using the Matlab routine and the results are for illustration 




Table 3.1: Calculation of the nonlinear terms for Section [3. 4. 21 



listed in Table O The resulting CMEs are 

(fi + a(d 2 yi + dl)) At + ( 7 o|^i| 2 + 2 7 i(|^ 2 | 2 + \A 3 \ 2 )) A, + l2 {A\ + ^Wi = 0, 

(fi + a(dl + d 2 y2 )) A 2 + ( l0 \A 2 \ 2 + 2 11 (\A 1 \ 2 + \A 3 \ 2 )) A 2 + l2 {A\ + Al)M = 0, (3.32) 

(fi + a{d 2 yi + d 2 J) A 3 + ( lQ \A 3 \ 2 + 2 11 (\A 1 \ 2 + \A 2 \ 2 )) A 3 + l2 (A 2 + A 2 )A 3 = 0, 
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where the following symmetries have been used: 70 := ^i. 1,1.1 = -^2.2,2.2 = ^3,3.3,3 due to ()3.250 with v = 7r/3; 
7i := ^1,2,2,1 = ^2,3,3,2 = -£1,2,1,2 = £2,3,2,3 due to 1)3.250 with v = n/3, and Q3.230 . Moreover, 71 = 7 2 ,i, 1,2 = 
^1.3,3.1 = -£3,1.1,3, where the second equality follows from (13.250 with v = tt/3 and the facts that few = 
r V3 (fc (3) " and u n {k^ - b^;x) = u n {k^;x) for all n e N. Finally 72 := I 2)2 ,i,i = £3,3,2,2 = = 
£2,2,3,3 due to (f3T25|) and ([3~24)) . and 72 = /2,2,1a = £l, 1,3,3 using (|3T25|) together with = r 7r/3 (fc (3) - 
and u n {kS z ^ — b^;x) = u n {k^;x) for all n e N. All the nonlinear coefficients are real: 70,71 £ R due to 
(|3.24|) and 72 G R since 72 = £2,2,1,1 = £1,1,2,2 by 1)3.240 and at the same time 72 = £2,2,1,1 = £1,1,2,2 by 1)3.260 . 
where we are using the facts that k^ — (k^\ —k^) T and that k^> = k^ does not hold. 

4. Numerical Examples of Gap Soliton Approximations. We compute here numerically localized 
solutions of the CMEs for the examples s 2 , S5 in Section 13. 4. II Then, using the leading order term in (|3. 1[) . we 
generate and plot an approximation of a gap soliton of the nonlinear Maxwell problem ()1.90 . In the evaluation 
of (|3.10 we position the photonic crystal so that the center of one of the annuli lies at the origin x = 0. 



4.1. Gap Soliton near the Edge s 2 . Figure Hj] plots in (a) the unique positive localized solution, the 



so called Townes soliton, of ()3.280 for the case w* = s 2 and in (b) the intensity I = \Ei\ 2 + \E 2 \ 2 + I-E3I 2 
of the leading order term in (|3. 10 . In Figure l4~2l we show the absolute value of the individual components 
Ex, E-2,E%. As the Townes soliton is radially symmetric, it was computed using the shooting method on ()3.280 
in polar coordinates. The fourth to fifth order explicit Runge-Kutta method ODE45 of Matlab was used in 
the shooting method. 

(a) A(y) W lE^IE/HE/ 
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Fig. 4.1: (a) CME solution, (b) intensity of the gap soliton approximation for the case u* 
See Section HH] 



4.2. Gap Soliton near the Edge S5. Here we restrict to solutions of (13.290 with the symmetry 

Ax = A A , A 2 = A 5 , A 3 = Aq, 

which reduces the problem to a system of three equations for Ax, A2, A3. To find a localized solution, we 
first replace fi by 0, and ax, a 2 , /?i, /3 2 by the average of these four numbers. Also the coefficients in each Mj, 
j G {1, 2, 3}, are replaced by their average. For this modified system the Townes soliton with Ax = A2 = A3 
is computed via the shooting method as in Section 14.11 Then a numerical homotopy in the coefficients is used 
to get a solution of Q3. 290 . The homotopy is applied to a fourth order centered finite difference discretization 
of (13.290 . Our homotopy always results in Ax = so that in the end we produce a solution of (13.290 with 
Ax = A4 = and A 2 — A 5 ^ 0, A3 — A 6 ^ 0. The two components A 2 ,A 3 are plotted in Figure |4"U1 together 
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(b) |E,1 




Fig. 4.2: Absolute value of the components Ei,E2,Ez of the gap soliton approximation for 
uj* = S2. See Section I4TT1 



with the intensity of the corresponding leading order term in (|3.1|> . In Figure FOl we plot the individual 
components of E in absolute value. 
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Fig. 4.3: (a) CME solutions A2, A3, (b) intensity of the gap soliton approximation for the case 
oj„ = S5. See Section \4. 21 



5. Conclusions. We have considered monochromatic out-of-plane gap solitons in Kerr nonlinear 2D 
photonic crystals as described by the full vector Maxwell system. Using a model of the nonlinear polariza- 
tion which does not produce higher harmonics, we arrive at a cubically nonlinear curl-curl problem for the 
fundamental harmonic. For gap solitons with frequencies in spectral gaps but in an asymptotic vicinity of a 
gap edge we assume a standard slowly varying envelope approximation based on the gap edge Bloch waves 
modulated by slowly varying envelopes of small amplitude. These envelopes are then shown to satisfy a system 
of coupled mode equations (CMEs) of the same structure as in the case of gap solitons of the 2D periodic 
nonlinear Schrodinger equation [13 ] I14 j . In particular the system generally involves mixed derivatives. Being 
a constant coefficient system depending only on the slow variables, the CMEs is a simple effective model for 
the near edge gap solitons. Similarly to [13] the derivation of CMEs needs to be carried out in Bloch variables 
due to the possible quasi-periodicity of gap edge Bloch waves. Symmetries among the coefficients of the CMEs 
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are determined using symmetries of the band structure and among the Bloch waves. 

We provide an example of a photonic crystal with a hexagonal periodicity lattice and a circular material 
structure in the periodicity cell. For this crystal three gaps are numerically observed (for u> > 0). CMEs are 
then derived for several gap edges including a case where a system of six CMEs arises. Numerical computations 
of localized solutions of these CMEs and of the corresponding gap soliton approximations are then performed. 
For the CME system with six components only solutions with four nonzero components were numerically 
constructed and it is unclear whether a solution with all six nonzero components exists. 

A rigorous justification of the CMEs, which states that for a certain class of CME solutions the full Maxwell 
system has gap soliton solutions which are indeed approximated by the slowly varying envelope asymptotic 
expansion, is expected to hold by similar arguments to those in [TJ1 [T31 [T3]. It will be the subject of future 
work. 
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